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I. INTRODUCTION 

The term Riemann problem usually refers to an initial value 
formulation for a hyperbolic set of partial differential equa- 
tions where initial data consist of two constant states separated 
by the discontinuity in the form of a plane surface. 

Although more than 150 years have passed since the pub- 
lication of Riemann' s original paper [1|, the significance of 
the Riemann problem in relativistic hydrodynamics has been 
truly recognized only recently (see [2| for an introduction). 

In this paper we are concerned with the general case of 
the Riemann problem in relativistic hydrodynamics, where the 
fluid is allowed to move with the velocity tangent to the initial 
discontinuity. In Newtonian hydrodynamics such solutions 
are not qualitatively different from those where the gas flows 
only in the direction normal to this surface. In the relativistic 
case, however, all velocities couple through Lorentz factors, 
and the resulting wave pattern of the solution depends on the 
velocities tangent to the initial discontinuity. Relativistic solu- 
tions of the Riemann problem with non-vanishing tangential 
velocities were first obtained numerically by Pons, Marti and 
Miiller in 2000 for the perfect gas equation of state [4|. Some 
solutions for other equations of state were also computed in 
IS El, and an analytic solution for the ultrarelativistic equa- 
tion of state was presented by Piftka and the author in [5 1. 

Solutions of this type are usually used in order to construct 
and test three dimensional hydrodynamical codes. However in 
2002 Aloy and Rezzolla discovered a hydrodynamical boost- 
ing mechanism accelerating gas in relativistic astrophysical 
jets to large Lorentz factors (3J- The analysis of this mecha- 
nism was based on solutions of the relativistic hydrodynami- 
cal Riemann problem with non-zero velocities tangent to the 
initial discontinuity. It is thus clear that the role of this class 
of Riemann problems is not restricted to numerical methods. 

This paper is devoted to three dimensional numerical stud- 
ies of corrugation instabilities occurring in solutions of the 
Riemann problem with non-zero tangent velocities both for 
ultrarelativistic and perfect gas equations of state. Such insta- 
bilities, mainly of the Kelvin-Helmholtz type, were suggested 
e.g. in [3]. Since solutions of the hydrodynamical Riemann 
problem usually consist of different elementary waves: rar- 
efaction and shock waves, a contact discontinuity and some 
constant states, the question of their stability is a little bit sub- 



tle. It makes sense to ask whether all parts of the solution 
become unstable or not. An instability of the region around a 
contact discontinuity does not have to imply that the surface of 
a shock wave present in the solution has to be unstable. Also, 
when performing studies of the stability of solutions, one has 
to decide on a certain notion of stability, and the final answer 
may depend on this choice. We will address these issues in 
Sec. V of this paper. 

We restrict ourselves to problems with moderately relativis- 
tic fluid velocities (v ~ 1/2), and concentrate on solutions 
characterized by different wave patterns, i.e., with two shock 
waves, two rarefaction waves or a shock and a rarefaction 
wave respectively. The question of stability of solutions with 
ultrarelativistic speeds is still open. 

There exist many numerical codes solving equations of rel- 
ativistic hydrodynamics for the perfect gas equation of state, 
or for general equations of state where the pressure can be ex- 
pressed in terms of the rest-mass density and the specific in- 
ternal energy. The case with ultrarelativistic equation of state 
is different in this respect; it basically requires a new formula- 
tion of the numerical scheme. Such a numerical code has been 
written for the purpose of this work and will also be described 
here. 

The order of this paper is as follows. In the next section we 
collect all basic equations necessary to explain our approach. 
Sec. Ill discusses the structure of exact solutions of the Rie- 
mann problem in relativistic hydrodynamics. In Sec. IV we 
describe new elements of numerical methods used in this pa- 
per. Sec. V contains numerical results on the stability of the 
Riemann problem in relativistic hydrodynamics. A short sum- 
mary can be found in Sec. VI. 



II. BASIC EQUATIONS 

Conservation of the energy and the momentum in relativis- 
tic hydrodynamics can be expressed as 

= o, (i) 

where T 1 ™ denotes the perfect fluid energy-momentum tensor 
T» v = (p+p)u»u v + p?f v . (2) 
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Here p is the pressure, p the energy density, u^ 1 denote compo- 
nents of the four-velocity of the fluid, and rf v are the compo- 
nents of the Minkowski metric. In the above formulae Greek 
indices range from to 3 and the convention with 

?f v = diag(-l,+l,+l,+l) 

is assumed. If the equation of state is in the barotropic form 
P - p(p)> m e above equations constitute a complete set of 
equations of hydrodynamics. The ultrarelativistic equation of 
state of the form p = c 2 s p, with c 2 being a constant (the local 
speed of sound) is a good example here. 

More generally one introduces the rest-mass density n, so 
that the continuity equation 

dft (nuf) = (3) 

is satisfied, and the equation of state is expressed in terms 
of the rest-mass density n and the specific internal energy e, 
defined as e = (p - n)/n. This is the case of the perfect gas 
equation of state: p = (y — l)«e, where y is a constant. In 
such a case the complete set of equations of hydrodynamics 
consists of four Eqs. ([TJ and Eq. ([3J. 

For the purpose of evolutionary problems, it is convenient 
to rewrite Eqs. ([T]) and Eq. ([3]) in the form, where the deriva- 
tives with respect to time and space are separated explicitly. 
We introduce the Lorentz factor W = u° and components of 
the three-velocity v' = u' jW. Eqs. ^ can be now written as 

d,U + fl,-P = 0, (4) 

where 



IE. SOLUTIONS OF THE RIEMANN PROBLEM 

Solutions of the Riemann problem for system I were ob- 
tained by Smoller and Temple for the case without veloci- 
ties tangent to the initial discontinuity in [9|. Solutions valid 
for arbitrary velocities (also tangent to the discontinuity) were 
presented by Pi§tka and the author in 0. 

Analogous solutions for system II were first calculated by 
Marti and Miiller for the case without tangential velocities [ 8 ] 
and promoted later to the case with non-zero tangential veloc- 
ities by Pons, Marti and Miiller in ||4). 

In order to describe these solutions briefly, we will intro- 
duce Cartesian coordinates x, y, z, and assume that the initial 
discontinuity, dividing two initial states L and R, is a plane 
given by x — 0. Here, traditionally, L refers to the "left" state, 
i.e., for x < 0; R refers to the"right" one, that is for x > 0. In 
this case solutions of the Riemann problem are self-similar: 
they depend on x and time t through a single variable £ = x/t, 
and they are independent of y and z. 

In all cases solutions can be constructed from self-similar 
elementary waves: a smooth part — the so-called rarefaction 
wave % and two possible discontinuities, i.e., a shock wave S 
and the contact discontinuity C. All three kinds of elementary 
waves can be separated by some constant intermediate states 
L„ and R t . Moreover, it can be shown that the initial state LR 
decays into one of the following possible wave patterns: 



V. 



U = (ip + p)W 2 -p,(p + p)W 
(p + p)W 2 v 2 ,(p + p)W 2 v 3 f 



(5) 



LR -> LS^L t CR,S^R, 
LR -> L'R+-L,CR»'R^R, 
LR -> LS^LtCRtft^R, 
LR -> LK-L,CR,S^R, 



F = ((p + p)W 2 v\ (p + p)W 2 v i v 1 + 6 n p, 

(p + p)W 2 v i v 2 + 6 i2 p, (p + p)W 2 v i v 3 + 6 a pf , (6) 

and (5' ; denotes the Kronecker's delta. 

The set of Eqs. ([T]l and Eq. Q can be also expressed in the 
form of Eq. To this end, it is customary to introduce the 
specific enthalpy h — (p + p)/n = 1 + e + p/n. Vectors U and 
F' must be now 5 -dimensional, namely 

U = [nhW 2 - p,nhW 2 v\nhW 2 v 2 ,nhW 2 v 3 ,nW) T (7) 

and 

F'' = (nhW 2 V, nhW 2 v'v l + 5 il p, 

nhW 2 v i v 2 + 6 i2 p, nhW 2 v'v 3 + S i3 p, nWv'f . (8) 

In this paper we will refer to Eqs. Q with U and F ! given 
by |5]l and |6]) and the ultra-relativistic equation of state as to 
system I. Eqs. ^ with U and F' defined as in (|7ji and ([8j and 
the perfect gas equation of state will be called system II. 



where the arrows refer to the direction in which the waves 
move with respect to each other (it is customary to refer to 
these waves as to left and right-moving). By the contact dis- 
continuity we understand a surface across which the pressure 
and the normal component of the velocity are continuous, 
while other quantities, like tangent components of the velocity 
or, in case of system II, the rest-mass density and the specific 
internal energy, can exhibit a jump. Accordingly, the pressure 
p and x components of the velocity v x have to be equal in both 
intermediate states L, and R*. Since it is possible to compute 
the relation between p and v x behind the left and right mov- 
ing waves (either a shock or a rarefaction wave) in terms of the 
state in front of the wave (L or R depending on the direction of 
propagation of the wave), we can compute the values of p and 
v x in both intermediate states from the condition that these val- 
ues have to be the same for both right and left moving waves. 
This calculation identifies the particular type of both waves. 
Next, it is possible to compute the speeds of propagation of 
discontinuities and of the front and the tail of the rarefaction 
wave (if present). At this stage computing all other details of 
the solution is straightforward. A precise description of this 
procedure can be found in JH |5] [8) . 
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IV. DESCRIPTION OF THE NUMERICAL CODES 

The code used to investigate the instabilities occurring in 
the solutions of the Riemann problem is based on a Godunov 
type, high resolution shock capturing scheme (for a textbook 
exposition see ifTOlfTTI '). The construction of the code is sim- 
ilar to that described in fl2l . and a preliminary version has 
been presented in |fT3l . 

The temporal evolution of conserved quantities is imple- 
mented according to a variant of method of lines. The space is 
discretised into cells (zones) labeled by indices i, j, k in Carte- 
sian directions x, y and z respectively. In the following x i7 yj 
and Zk will denote coordinates of the center of the zone labeled 
by i, j and k. Dimensions of the zone will be given by Axi, 
Ay, and Az&. Positions of the interfaces between zones will 
be denoted in the usual fashion, where symbol x I+ i/2 refers to 
the x coordinate of the interface between zones labeled by i 
and i + 1, and positions of the interfaces in directions y and z 
are denoted in the analogous way. The time derivative of cell- 
averaged values of U is computed according to the following 
formula 



dV 



U,k 



dt 



~ Ax/ ft' +l / 2 'j' k *"i-V2,j,k) 



—L(p v -P ) 

Ay . \ r iJ+U2,k r iJ-l/2,k) 

~~£z~ k fttjMW ~ ^lj,k- 1/2) ' 



where F r , , . , , F A ' , „ . , , F* , W. . , „ , , F z . , , , „ and 
F z . , , n denote numerical fluxes defined at the interfaces be- 
tween adjacent cells. Values of U, ^ are advanced in time us- 
ing the standard fourth or second order Runge-Kutta method, 
i.e., either as 

u "ji = u u* + l At ^ + 2k 2 + 2k 3 + k 4 ) , 



where 



ki 
k 2 
k 3 
k 4 



dVij,k 
dt 



dt 

dV u ,k 
dt 



d\Ji jk 1 „ \ 



or according to 



dt 



Afk 2 . 



Here the upper index n numbers subsequent time steps of size 
At. 

Clearly, the key problem is to compute the values of nu- 
merical fluxes F. In [ 13 1 we have presented preliminary tests 
obtained by using relativistic Harten, van Leer, Lax, Ein- 
feldt (HLLE) formulae for system I (for the description of the 



HLLE solver see |[T4l4T6ll ). Here we go a little bit further and 
employ a modified version of a flux formula introduced origi- 
nally by Donat and Marquina 1 17 1 and used in |[T2l . 

Suppose we want to compute a numerical flux at the inter- 
face between two Riemann states Ul and Ur. Let A p , l p and 
r p be the eigenvalues, left eigenvectors and right eigenvectors 
of the Jacobian <9F' /3U respectively. Moreover, subscripts L 
and R will refer to the values obtained for the left and right 
states Ul and XJr. Numerical fluxes F' are computed as 

P = l - (F ! (U L ) + F ! (U«) 

- ^ max L , R |^ p | (q. pJ i ■ U«)r p>R - (l p>L ■ U t )r p , £ ) i . 
p ) 

The values of A p , l p and r p should be computed analyti- 
cally. They are given in [18| for system II, and have to be 
computed separately for system I. It can be noted, however, 
that the precise knowledge of all these terms is not necessary. 
The only required quantities are the eigenvalues A p and all 
vectors (l p ■ \J)r p . 

In general, for a barotropic fluid with p = p(p) and c 2 s - 
dp I dp the eigenvalues of dF' /d\J are 

A = v'', 



A + = 



v'd - c]) ± c s V(l - v k v k ) (1 - v*v*c5 - (v0 2 (l - c?)) 



1 - \'kV k C 2 



The eigenvalue A$ is twofold degenerate. The expression for 
A ± can be also written as 



A ± = 



v> ±A 
1 + v'A 



with A' 2 = 1 + W 2 (l - (v ! ) 2 ) (1 - c 2 )/c 2 s . If V is the only 
component of the velocity then A — c s . 

The remaining formulae will be given for dF x /d\J. Expres- 
sions for SF V /(9U and dF : /d\J can be easily obtained from the 
symmetry. (There is also no reason to treat Jacobians dF' /d\J 
separately in the code. It is enough to implement formulae for 
F A ; other fluxes can be computed by swapping the order of ve- 
locities and conserved momenta.) Let us introduce the square 
of the tangential velocity vf = v 2 + v 2 , and define a bunch of 
auxiliary quantities: 

® ± =(? s p(\-v 2 j(^-J x -\+2A*v x ), 

E± = (1 - *±v x ) (v 2 (l - v 2 ) + 2v) - v 2 (l - v 2 )) 

-2v 2 (l-v 2 ), 

Q^^a + ^v^-a+c 2 ^, 

" ± (® ± + p(c 2 s Z ± + ( 1 - v k v k )(A T - v x )v x j) 



In these terms 

2 



pW 2 



Yjdoj ■ U)r 0J = f—; (2v,v 2 , v,,(l - v 2 . + v 2 ), 
v.(l -v 2 + v 2 ),2v 2 ) r 



7=1 
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and 

(1 ± • U)r ± = 

((A±~ v x )v x + c] (l - vf - A ± v x (2 -v 2 x + vf)) 
{ ^(l + ^-^-fl+av, 
Vx (c]{^-^)-\) + A ± {\-c]v k v k W 
V) " VZ ' i ± (l+c?(v?-v?))-(l+ C ?)v, , 

Systems I and II differ also in the numerical procedure used 
to recover primitive hydrodynamical quantities like n, p, v' 
from the conserved ones. For the perfect gas equation of state 
such a recovery is performed by means of a Newton-Raphson 
scheme. For the ultrarelativistic equation of state primitive 
values can obtained from U analytically. 

The last, key ingredient of the code is the reconstruction 
procedure used to obtain the states Ul and Us based on the 
values of U from the neighboring zones. For the results pre- 
sented here we have used the Convex Essentially Non Oscilla- 
tory (CENO) reconstruction procedure presented in |[T9l and 
applied to the relativistic hydrodynamics by [20, l2TI . 

The general rules of the method can be found in [ 19 1. Here 
we will only give final formulae. The CENO reconstruction 
is applied in each dimension to each component u of the con- 
served vector U separately. Consider the x direction and a 
zone Xi-ip < x < x i+l /2- Let us define 



St 



(Ax;) z u i+l - ((Axfy - (Ax;) z ) Ui - {AxfYu, 
AxfAxJ(Ax? + Ax~) 

AxJUi+i — (Ax? — Axj)ui — Ax^Uj-i 



Si = 2 



Ax? Ax. (Ax + + Ax. ) 



and 



Si = mm | : — ; ,S i 



Ax? 



Ax. 



where Ax? = x-, + \ - x,-, Ax i = x-, - Xi-\, and mm denotes the 
minmod function, i.e., 

|min{xi , . . . , x n }, if x, > for all i, 
max{xi ,...,x n }, if x, < for all i, 
0, otherwise. 

Next, we introduce a linear function 

Li(x) = Ui + Si(x - xi) 

and three quadratic polynomials 



Qf } = Ui+k +S i+k (x - Xi) + ^S i+k (x - XiY, 

fork - 0, +1. Following |19| we will now choose quadratic 
polynomials gP that are closest to L(x) at x,_i/2 and Xi+i/2 
respectively. More precisely, in order to obtain the left state 
u f+i/2' we com P ute differences d {li} = Qp (Xi+1/2) - ^>(JCj+i/2). 
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TABLE I. Initial data for sample problems that were evolved numer- 
ically. Problems (a)-(c) correspond to equations of system I. Prob- 
lems (d)-(f) apply to system II. The resulting wave pattern is indi- 
cated as "type". Symbols SS, KR and RS refer to configurations 
with two shock waves, two rarefaction waves and a combination of a 
shock and a rarefaction wave respectively. 

System I 



Problem 


PL 




V 

V L 


Pr 




V 

V R 


Type 


(a) 


0.5 


0.2 


0.2 


0.5 


-0.2 


-0.2 


SS 


(b) 


0.5 


-0.2 


0.2 


0.5 


0.2 


-0.2 


KR 


(c) 


0.5 





0.2 


1.0 





-0.2 


US 


System II 


Problem 


n L £l 


vZ 


V 

V L 


»« £r 


v "r 


v l 


Type 


(d) 


1.0 0.5 


0.2 


0.2 


1.0 0.5 


-0.2 


-0.2 


SS 


(c) 


1.0 0.5 


-0.2 


0.2 


1.0 0.5 


0.2 


-0.2 


KR 


(0 


1.0 0.5 





0.2 


1.0 1.0 





-0.2 


US 



In the case where all S k) have the same sign we choose the 
polynomial Qi = Q.f^ for which |c/ (i:) | is the smallest. Other- 



wise we select the linear function Q, = L, . The left state u L 



is now obtained as uf +1 ^ 2 



1+1/2 



Qi(xi+\/2)- An analogous selection 
procedure is applied at the interface x,_i/2 in order to compute 

«f_l /2 = Qi(Xi-l/l)- 

The code was tested against different exact solutions of the 
Riemann problem. The results were slightly better than those 
presented in Ifl3ll . We have also performed convergence tests 
for three dimensional runs with satisfactory results, i.e., the 
solution computed on a coarser grid resembled that computed 
with a better resolution. 



V. CORRUGATION INSTABILITY PROBLEM 

In this paper we investigate the stability of the relativistic 
Riemann problem with non-zero tangent velocities. We will, 
however, restrict ourselves to corrugation instabilities only — 
perturbations of the initial data will be applied to the shape of 
the surface dividing two initial states. 

Stability of simple waves in the solutions will be first under- 
stood according to a very simple notion proposed by Anile and 
Russo in Il22ll . Their idea can be summarized as follows: con- 
sider a corrugated shock wave propagating into some medium. 
If a convex part of the perturbed shock wave moves with a 
speed that is larger as compared to the speed of an unper- 
turbed shock wave, and reciprocally, if a concave part propa- 
gates slower than an unperturbed shock wave, then the pertur- 
bations would tend to amplify, and we would say that the wave 
is unstable. A converse situation would lead to a smoothing 
of the shock wave and a decrease of the size of corrugations. 
Such behavior is called stable. 

This clearly geometrical notion of stability can be accom- 
panied with some global quantitative data. In order to measure 
the size of evolving perturbations we have compared each per- 
turbed solution with a corresponding unperturbed one (i.e., a 
solution for which the initial discontinuity is a plane surface 
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FIG. 1. Three dimensional distribution of conserved energy density e 
for problem (a). Subsequent snapshots correspond to evolution times 
t = 0.5,1.0, 1.5,2,2.5. 



and both initial states are the same as in the original solution). 
We have decided to look at the perturbation in the conserved 




FIG. 2. Three dimensional distribution of conserved energy density e 
for problem (b). Subsequent snapshots correspond to evolution times 
t = 0.4,0.8, 1.2, 1.6,2.0. 



energy e — (p + p)W 2 - p. Technically, we introduce the 




for problem (c). Subsequent snapshots correspond to evolution times 

t = 0.0, 0.6, 1.2, 1.8, 2.4. FIG. 4. Three dimensional distribution of conserved energy density e 

for problem (d). Subsequent snapshots correspond to evolution times 
t = 0.9, 1.8, 2.7, 3.6, 4.5. 
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FIG. 5. Three dimensional distribution of conserved energy density e 
for problem (e). Subsequent snapshots correspond to evolution times 
t = 0.6, 1.2, 1.8,2.4,3.0. 



FIG. 6. Three dimensional distribution of conserved energy density e 
for problem (f). Subsequent snapshots correspond to evolution times 
t = 0.0,0.9, 1.8,2.7,3.6. 
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FIG. 7. Evolution of the rest-mass density n for problem (f). The 
plots show two dimensional z = const cross sections through the 
grid. Subsequent snapshots correspond to evolution times t = 
0.06,0.9,1.8,2.7,3-6. 



FIG. 8. Evolution of the v v component of the velocity for problem 
(f). The plots show two dimensional z = const cross sections through 
the grid. Subsequent snapshots correspond to evolution times t = 
0.06,0.9,1.8,2.7,3.6. 



following three L norms: 

||^ ^unperturbed) | ^ i — ^ ^ ^ 



- unperturbed h 



'Unperturbedj | t2 



cPi 



\& ^unperturbedl' 



j|^ ^unperturbedjl^ — SUp|£ ^unperturbed I ! (9) 

and we will be interested in their evolution in time. 

The Riemann problem is evolved numerically on a 3- 
dimensional Cartesian grid of 1800 x 600 x 300 zones span- 
ning a cuboid region x e [-1.5, 1.5], y € [0, 1], z e [0,0.5]. 



The initial discontinuity is located at the center of the grid, 
and, modulo perturbations, it is directed perpendicular to the 
x axis. 

On the boundaries in directions y and z, that is on the faces 
perpendicular to the initial discontinuity, we assume periodic 
boundary conditions. On the remaining two faces, that is in x 
direction, outflow boundaries are implemented. 

Perturbations of the surface of initial discontinuity are ap- 
plied according to the following formula 

0, r > R, 
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FIG. 9. The growth of the perturbation in the energy e computed with 
respect to L 1 and L 2 norms (Eqs. for problem (a). 



FIG. 11. The growth of the perturbation in the energy e computed 
with respect to L 1 and L 2 norms (Eqs. {9}) for problem (c). 
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FIG. 10. The growth of the perturbation in the energy e computed 
with respect to 1} and L 2 norms (Eqs. |9]l) for problem (b). 



FIG. 12. The growth of the perturbation in the energy e computed 
with respect to L 1 and L 2 norms (Eqs. {9}) for problem (d). 



x(y,z) 



(x - xif + (y- yi) : 



Here A,, x,, y,- are random amplitudes, radii, and coordi- 
nates of the center for each particular perturbation. A little 
care is required in order to make those perturbation compati- 
ble with the assumption of periodicity at the boundaries. 

We decided to focus on a set of six solutions (a)-(f) with 
moderately relativistic fluid velocities (v ~ 1/2), but exhaust- 
ing all possible different wave patterns for systems I and II. 
Values characterizing initial states for these problems are col- 
lected in Table 1. System I was evolved with c 2 s = 1/3. For 
system II we have chosen y — 4/3. This corresponds to equa- 
tions of state p = p/3 and p = ne/3 respectively. 

Graphs showing the evolution of the initial data for all prob- 
lems (a)-(f) are presented on Figs. [TJj6] In each case a series 
of three dimensional perspective plots corresponding to sub- 
sequent moments in time is shown. The surfaces on the plots 
are the surfaces of constant energy density e (isopycnic sur- 
faces). Also, on the back face of the grid the boundary values 
of the energy density e are color-coded. All pictures are ori- 
ented in such a way that the x axis is pointing slightly in the 
direction of the viewer. 



0.08 
0.07 
0.06 
0.05 
0.04 
0.03 
0.02 
0.01 




FIG. 13. The growth of the perturbation in the energy e computed 
with respect to L 1 and L 2 norms (Eqs. |9}) for problem (e). 



The common feature that can be observed on all plots is the 
existence of an unstable turbulent region around the contact 
discontinuity with instabilities in the form of twisted tubes of 
considerably lower density. Another common feature is that 
in all cases both shock waves and rarefaction waves behave in 
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FIG. 14. The growth of the perturbation in the energy e computed 
with respect to L 1 and Lr norms (Eqs. {9]l) for problem (f). 
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FIG. 15. The growth of the perturbation in the energy e computed 
with respect to L°° norm (Eqs. |9j) for problems (a), (b) and (c). 



The instability occurring around the contact discontinu- 
ity can be identified as a relativistic version of the Kelvin- 
Helmholtz instability. In all cases (a)-(f) the gas in the in- 
termediate states L* and R* flows with some nonzero tangent 
velocity v y , different at both sides of the contact discontinu- 
ity. We have checked that in the absence of the tangent ve- 
locity no instability around the contact discontinuity devel- 
ops. The Kelvin-Helmholtz character of the instability can be 
seen on Fig. [7] presenting once again the evolution of prob- 
lem (f). Here, we have decided to show the distribution of 
the rest-mass density n on two dimensional z = const cross 
sections through the grid. In this case the two intermediate 
states L, and R* have different densities and we can observe a 
standard Kelvin-Helmholtz wave picture. The tubes of lower 
density can be identified with regions trapped under crests of 
the waves. 

In the same fashion Fig. [8] shows the two dimensional plots 
of the v y component of the velocity for problem (f), i.e., for 
the same data as those presented on Fig. [7] 

Figs. [9 16 show the time evolution of norms (|9]l for so- 
lutions (a)-(f). The evolution of L 1 and L 2 norms for those 
problems is depicted on Figs. [9 14 In both norms there is a 
short phase of a rapid growth of the perturbation followed by 
either a decrease of its size (L 2 norm) or a very slow growth 
(L 1 norm). For completeness L°° norms for problems (a)-(f) 
are shown on Figs. 15 and 16 Also this norm increases at 



early stages of evolution, and starts to oscillate around some 
constant value for later times. Judging from those results, 
one can say that the fast exponential growth of instabilities, 
usually obtained in a linearized calculation for the non rela- 
tivistic Kelvin-Helmholtz instability, is only limited to a very 
short early phase of evolution. Later, nonlinear effects become 
dominant, and the growth of perturbations is stalled. 



VI. SUMMARY 
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FIG. 16. The growth of the perturbation in the energy e computed 
with respect to L°° norm (Eqs. |9|) for problems (d), (e) and (f). 



a stable way — their shape gets smoother and smoother with 
time. One can notice, however, that the behavior of problems 
with the perfect gas equation of state is more turbulent than 
that of problems with the ultrarelativistic equation of state. 



We have performed a series of three dimensional numerical 
studies of the corrugation stability of the Riemann problem in 
relativistic hydrodynamics with non-zero velocities tangent to 
the surface of the initial discontinuity. We specialized to two 
equations of state: the ultrarelativistic one and that of perfect 
gas. In both cases a modern high resolution shock capturing 
Godunov type numerical scheme was employed. The con- 
served hydrodynamical quantities were evolved in time within 
the framework of method of lines, using a standard fourth or- 
der Runge-Kutta algorithm. Numerical fluxes were based on 
the work of Donat and Marquina IfTTl . The required spectral 
decomposition of the Jacobians appearing in the equations of 
relativistic hydrodynamics was taken from fTHl for the per- 
fect gas equation of state; it was computed separately for the 
general barotropic case (ultrarelativistic equation of state, in 
particular), and the results are presented in this paper. The 
CENO reconstruction procedure was adapted from lfT9Tl2T1l . 

We have focused on mildly relativistic solutions resulting 
with different possible wave patterns, i.e., solutions with two 
shock waves, two rarefaction waves, or a combination of a 
shock wave and a rarefaction wave. In all cases both rarefac- 
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tion and shock waves behaved in a stable way in the sense 
that their shapes were becoming flattened with time. Kelvin- 
Helmholtz type instabilities were observed to develop only 
around the contact discontinuity, forming a turbulent region 
with characteristic rarefaction "tubes". This behavior is es- 
sentially the same for the perfect gas equation of state and for 
the ultrarelativistic one. Slightly more turbulent instabilities 
develop in the case of perfect gas equation of state, what can 
be observed both on the three dimensional plots of solutions 
and in the evolution of norms measuring the size of the per- 
turbation. The behavior of perturbed solutions with ultrarel- 
ativistic velocities tangent to the initial discontinuity remains 
an open question. 

We believe that these results can be helpful in understand- 
ing more complex phenomena occurring, for instance, during 
propagation of relativistic astrophysical jets. The evolution of 
the so-called turbulent cocoon [23 1 could be compared with 
our results concerning the Riemann problem alone. Since the 
behavior of turbulent flows depends on dimensionality, and in 



particular it is different in two and three spatial dimensions, a 
fair comparison would require three dimensional simulations 
of jets (see [24 1 for an example). 

Solutions of the relativistic Riemann problems are also ap- 
plied in the modelling of heavy ion collisions during the so- 
called hydrodynamical phase (cf. ||251 and references therein). 
Our results should be of interest here. 
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